function [BB, BB_deriv_1, BB_deriv_2] = construct_splines_twodim(x_interest_1, x_interest_2, xi_1, xi_2, n1, n2, k1, k2, p_low1, p_low2, p_high1,p_high2, down)

[B_1, B_deriv_1]=construct_splines(x_interest_1, xi_1, n1, k1, p_low1, p_high1, down);
[B_2, B_deriv_2]=construct_splines(x_interest_2, xi_2, n2, k2, p_low2, p_high2, down);

BB=[];
BB_deriv_1=[];
BB_deriv_2=[];

for ii=1:(n2+k2),

    B_add=B_1.*(B_2(:,ii)*ones(1,n1+k1));
    BB=[BB B_add];    

    B_deriv_add1=B_deriv_1 .* (B_2(:,ii)*ones(1,n1+k1));
    BB_deriv_1=[BB_deriv_1 B_deriv_add1];
    
    B_deriv_add2=B_1 .* (B_deriv_2(:,ii)*ones(1,n1+k1));
    BB_deriv_2=[BB_deriv_2 B_deriv_add2];
    
end;

    %
